Early assembly of the most massive galaxies 
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The current consensus is that galaxies begin as small density fluctuations in the early 
Universe and grow by in situ star formation and hierarchical merging ^ Stars begin 
to form relatively quickly in sub-galactic sized building blocks called haloes which 
are subsequently assembled into galaxies. However, exactly when this assembly takes 
place is a matter of some debate^' ^. Here we report that the stellar masses of brightest 
cluster galaxies, which are the most luminous objects emitting stellar light, some 9 bil- 
lion years ago are not significantly different from their stellar masses today. Bright- 
est cluster galaxies are almost fully assembled 4 — 5 Gyrs after the Big Bang, having 
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grown to more than 90% of their final stellar mass by this time. Our data conflict 
with the most recent galaxy formation models '''^ based on the largest simulations 
of dark matter halo development ^. These models predict protracted formation of 
brightest cluster galaxies over a Hubble time, with only 22% of the stellar mass as- 
sembled at the epoch probed by our sample. Our findings suggest a new picture in 
which brightest cluster galaxies experience an early period of rapid growth rather 
than prolonged hierarchical assembly. 

Brightest cluster galaxies (BCGs) are located at the centres of galaxy clusters. They 
constitute a separate population from bright elliptical galaxies ^ and both their homogene- 
ity and extreme luminosity have motivated their use as standard candles for cosmology 
Our investigation focuses on BCGs in the most distant X-ray emitting galaxy clusters at 
redshifts z = 1.2 — 1.5, where (l + z) is the expansion factor of the Universe relative to the 
present. It has been shown that X-ray cluster selection is currently the optimum strategy 
for an unbiased investigation of BCG evolution Properties of our BCGs and their host 
clusters are listed in Table 1. All five clusters were discovered serendipitously in X-rays 
and they are the most distant clusters discovered in their respective X-ray surveys 
The cluster J2215 was discovered as part of the XMM Cluster Survey (XCS and has 
the highest redshift of any spectroscopically confirmed cluster 
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The stellar mass of a BCG depends upon the hierarchical build up of its host dark 
matter halo and its stellar evolution history, along with the baryonic physics of the galaxy. 
We base our study of BCGs on photometry in the infrared wavebands J (1.26 /xm) and 
Kg (2.14 jim). Infrared imaging is essential at these large redshifts to compensate for the 
redshifting of the early-type galaxy spectra. Also, these wavebands are less sensitive than 
optical light to the presence of young stars and are a more accurate tracer of the underlying 
old stellar population and, hence, of the stellar mass of the systems. Fig. 1 shows an 
infrared image of the cluster J2235 from our sample (see also Supplementary Fig. 1). 

We start by examining the ages of the stars themselves in these galaxies using the 
run of J — Kg colour evolution with redshift as shown in Fig. 2. For BCGs at the redshift 
of our sample the J — Kg colour predictions for the models separate clearly. For the 
comparison sample at lower redshift we use X-ray selected clusters which are well 
matched in mass to our own cluster sample. There is a remarkable agreement between the 
data and the hybrid model (see Fig. 2 legend), with all five BCGs lying within 0.05 mags 
of their predicted colour, indicating a consistent epoch of formation for the majority of the 
constituent stars in all systems between redshifts zj = 3 — 5, some 2-3 Gyr after the Big 
Bang. 

Turning our attention to the mass assembly of BCGs implied by our data, in Fig. 3 
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(see also Supplementary Table 1 and Supplementary Fig. 2) we show the estimates of 
stellar mass for our distant BCGs normalised to the average mass of the comparison sample 
atz < 0.04, which is 8.99 (±0.82) x 10^^ Mq (s.e.m.). Using a Tukey's biweight location 
estimator for robustness, for our five objects located at 2; = 1.22 — 1.46 we find an average 
stellar mass of 8.86 (±1.73) x 10" Mq (s.e.m.). The ratio of these estimates is 0.99 ±0.21 
(s.e.m.), indicating that on average the masses of the high redshift BCGs are consistent 
with local counterparts. 

To compare with theory we use the haloes from the Millennium Simulation^ 



(http://www.mpa-garching.mpg.de/millennium) matched to the total mass of our clusters. 



estimated from their X-ray luminosity (see Supplementary Information). The mass range 
of our five clusters (Table 1) has excellent overlap with the combined 2; = 1.08 and z = 1.5 
halo samples ^ (Supplementary Fig. 3). The predicted hierarchical mass build up of BCGs 
in these 250 haloes is also shown in Fig. 3. The corresponding mass of the simulated 
BCGs has grown to an average of only 1.92 (±0.38) x 10^^ M© (s.d.) by this time, some 
22% of the observed value. The data are inconsistent with the prediction at the level of Aa 
(one-tailed P = 0.008, degrees of freedom d.f. = 4; based on a Student's t distribution 
appropriate for small samples). 

To check the stability of the BCG assembly predictions we selected massive haloes 
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from the independent Durham semi-analytic model ^ which also uses the Millennium 
Simulation^ but incorporates a different treatment of the baryon physics close to active 
galactic nuclei, partly in order to better reproduce the abundance of massive elliptical 
galaxies at high redshift. Using the same selection limits we find that the BCG mass 
fractions compared to the present day are 0.22;tao9 

at 2 = 1 and 0.1712 ;q7 at z — 1.5, 
indicating good agreement between the two semi-analytical models. 

It is well known that the estimates of stellar mass from photometry even for early- 
type galaxies such as BCGs depend on the underlying stellar evolution model used. To 
investigate this sensitivity we have applied three independent stellar population synthesis 
codes to early-type galaxies at the mean redshift of our sample (z — 1.3) using a range 
of model parameters (see Supplementary Table 1). These results show that the Kg band 
stellar mass estimates remain significantly discrepant from the semi-analytic predictions 
(one-tailed P < 0.02, d.f . = 4) for the vast majority of parameters considered across the 
three models, reaching a value for one-tailed Pof > 0.05 in one of the three only if the 
stellar formation epoch Zf is less than 2.5 together with a stellar metallicity less than the 
solar value. This situation is incompatible with observations of BCGs and massive early- 
type galaxies in general (see Supplementary Information). We conclude that there remains 
a significant discrepancy between the recent semi- analytic models of galaxy formation 
coupled to the largest N-body simulations and the stellar masses of BCGs at the centres of 
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the most massive clusters. 

In comparison to recent studies this work significantly extends the redshift base- 
line over which BCG evolution has been investigated to z = 1.5, equivalent to a look-back 
time ~ 65% the age of the Universe. Although the first glimpse of the z > 1 BCG popu- 
lation reveals galaxies with a range of stellar masses, there is on average considerably less 
stellar mass evolution than expected, with the bulk (> 90%) of the stellar mass already 
in place by z ~ 1.5, corresponding to only about 4-5 Gyrs after the Big Bang; the cur- 
rent models predict a considerably longer timescale of about 1 1 Gyr for the same growth, 
reaching 90% at 2; ~ 0.2. 

Despite this there is evidence that merging is still underway in our high redshift 
sample. The BCG in J0849 at ^ = 1.26 has a nearby companion (projected separation of 
about 6 kpc) with which it is likely to undergo dissipationless merging in the future Of 
the other clusters in our sample, the BCG and its neighbour (projected separation of about 
15 kpc) in J 1252 are also possible merger candidates. Assuming that mergers take place 
in both these cases, the fraction of BCG stellar mass already assembled (based on the Kg 
fluxes of the main components) is ~ 84% and ~ 60% for J0849 and J 1252 respectively, 
supporting the contention that most of the growth has actually already taken place in these 
two BCGs. 
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The timescale for the mass assemblage is similar to the age of the component stars 
(2 — 3 Gyrs), a situation that appears to resembles classical monolithic collapse ^^'^^ rather 
than hierarchical formation. To form a galaxy of stellar mass 10^^ Mq over 4 Gyrs re- 
quires a mass deposition rate of about 250 M© yr~^ and an efficient mechanism to feed 
the gas into the inner regions of the halo where it can form stars. Unfortunately the merg- 
ing process becomes inefficient for massive galaxies because merger induced shocks lead 
to heating as opposed to radiative cooling of the gas ^. One recent suggestion is that 
the early assembly of massive galaxies at 2; > 2 is driven by narrow streams of dense 
cold gas which penetrate the shock-heated region greatly increasing the efficiency of the 
gas deposition and associated star formation. Thus, the fraction of time that young BCGs 
spend undergoing a major merger event could be < 10%, with the stellar mass assem- 
bly dominated by this 'stream-fed' process Alternatively, a deficiency may lie in the 
semi-analytic treatment of the physical processes in the densest environments during early 
hierarchical assembly; this contention is supported by the fact that current predictions are 
moderately consistent with observations of the evolution of luminous red galaxies 
whereas our results, which focus on the most massive subset of this population, the BCGs, 
differ much more from the model predictions. 

In a wider context the hierarchical simulations and their semi-analytic prescriptions 
have arguably provided an excellent way of generating mock catalogues of galaxies to 
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compare with real data, but our results show that they do not account for the assemblage 
history of all galaxies. Larger simulations may provide a better statistical probe of both 
the merging history of the largest haloes and cluster-mass trends. If BCGs collapsed and 
formed at high redshift in a single burst of intense star formation then they may well be 
dusty and in sufficient numbers to be detectable with the coming generation of submil- 
limetre surveys, which will cover areas large enough to detect objects as rare as BCGs. 
The ongoing XCS survey will find many more high redshift clusters and we anticipate that 
our results will stimulate independent studies of BCGs as new clusters are found in the 
redshift 'desert' beyond z — 1.5 from infrared and X-ray based surveys such as eRosita. 



9 



Supplementary Information is linked to the online version of the paper at www.nature.com/nature. 

Acknowledgements This work is based in part on data collected at the Subaru Telescope, which 
is operated by the National Astronomical Observatory of Japan and the XMM-Newton, an ESA 
science mission funded by contributions from ESA member states and from NASA. We acknowl- 
edge financial support from Liverpool John Moores University and the STFC. M.H. acknowledges 
support from the South African National Research Foundation. IRAF is distributed by the Na- 
tional Optical Astronomy Observatories, which are operated by the Association of Universities for 
Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. 
We thank Gabriella De Lucia for making simulation results available to us in tabular form, Ichi 
Tanaka for developing the MCSRED package used to reduce the MOIRCS data, Maurizio Salaris 
for discussions on stellar population synthesis models and Ben Maughan for discussions on cluster 
masses. 

Author Contributions C.A.C. provided the scientific leadership, helped design the experiment, 
wrote the paper and led the interpretation. J.P.S. performed the photometry and data analysis and 
made major contributions to the interpretation. M.H. wrote the Subaru proposal, carried out the data 
reduction and photometric cahbration, contributed to the analysis and interpretation and provided 
detailed comments on the manuscript. S.T.K. independently checked the cluster mass calculations. 
S.A.S. provided useful discussions on the data and comments on the manuscript. The remaining 
authors make up the team of the wider XCS project which led to the discovery of J2215. R.G.M, 



10 



R.C.N, and A.K.R. made useful comments on the text. 

Competing Interests The authors declare that they have no competing financial interests. 

Correspondence Correspondence and requests for materials should be addressed to 
C.A.C. (email: cac@astro.livjm.ac.uk). 



11 



Figure 1. Infrared image of the cluster J2235. An infrared image of the cluster J2235 at a 
redshift z = 1.39. Data were taken using the 8.2-m Subaru telescope. The image is combined from 
separate J and exposures and shows the 1.5' x 1.5' region surrounding the cluster centre. At 
this redshift 1.5' corresponds approximately to 0.75 Mpc. The green overlaid contours show the 
smoothed X-ray emission taken from the XMM -Newton XCS pipeline, smoothed with a Gaussian 
kernel. The X-ray peak coincides with the cluster centre and the position of the BCG. For a full 
description of the observations and data reduction see Supplementary Information. 
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Figure 2. The stellar evolution of BCGs with redshift. The J — Kg colour evolution for our 
five high redshift BCGs (red) and 72 BCGs from the comparison sample (black) which have 
host cluster masses in the same range as our high redshift clusters and have available J and Kg 
photometry. The errors (s.d.) reported for the comparison sample and our data are ~ 0.1 
mag and ~ 0.02 mag respectively and are shown in the figure. This plot includes simple stellar 
population models incorporating: no stellar evolution (solid); passive evolution with formation 
epoch Zf = 5 (dashed); passive evolution with formation epoch Zf = 2 (dotted); a hybrid model 
with an exponentially decaying star formation rate in which 50% of the BCG stellar content is 
formed by = 5 and 80% by Zf = 3 (dot-dashed), which is appropriate to the star formation 
history predicted by the semi-analytic model The Zf = 2 and Zf = 5 stellar models are calculated 
assuming solar metalicity and a Salpeter initial mass function (IMF) while the hybrid model was 
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calculated with a Chabrier IMF ^. The implied epoch of formation 2;/ = 3 — 5 (2 — 3 Gyrs after 
the Big Bang) agrees well with other estimates of stellar ages determined for BCGs and early- 
type galaxies in clusters (see Supplementary Information). Throughout our analysis we assume a 
concordance cosmology of Clm = 0.3, Cl\ = 0.7, and Hq = 70 km Mpc~^, where ^Ia is the 
energy density associated with a cosmological constant. See the Supplementary Information for 
details of data reduction. 
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Figure 3. The mass evolution of BCGs with redshift. The BCG mass estimates of our sam- 
ple normalised to local galaxies at z = 0.04. The red cross is the estimated biweight location 
(8.86 X 10^^ Mq) and scale (3.87 x 10^^ M©) of the sample. We calibrate the stellar masses 
by comparing the rest-frame absolute Kg magnitudes with the predicted magnitudes and corre- 
sponding stellar masses from the semi-analytic models ^. This involves correcting the observed 
Kg for: cosmological dimming; sampling different spectral regions of the galaxies resulting from 
the redshift (fc-correction); and stellar evolution. The last two corrections are carried out using 
synthesized stellar spectra for early-type galaxies (appropriate to BCGs) from the hybrid stellar 
population model shown in Fig. 2. The A;-correction is well understood over the wavelength range 
appropriate to our sample (0.9 — 2.2 ^um) introducing an uncertainty of about 10% in the rest-frame 
absolute Kg magnitude estimates. The biweight scale provides a realistic estimate of the intrinsic 
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error (s.d.) in the average mass using the hybrid model, however the total uncertainty in the in- 
ferred BCG mass is larger as it depends on the stellar evolution model used (see Supplementary 
Information). The grey diamonds show the individual BCG mass predictions in 125 simulated 
clusters at each of six redshifts (0.0, 0.2, 0.5, 0.75, 1.08, 1.5) above corresponding selection masses 
(4.7, 3.5, 2.8, 2.4, 1.5, 1.0) in units 10^^ M©. The black filled circles show the average value at 
each redshift (all errors are s.d.). The predictions are based on semi-analytic models of galaxy evo- 
lution. These use large N-body simulations such as the Millennium Simulation \ which models 
the development of 2160^ cold dark matter particles within a box over 2 billion hght years on each 
side. The semi-analytic techniques use the merger trees from the simulations and graft on analyt- 
ical approximations to account for the comphcated physics of the baryons in a range of ongoing 
processes associated with galaxy formation, such as: cooling, star formation, supernova outbursts 
and the growth of black holes in active galactic nuclei. 
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Table 1: The properties of the host clusters and their BCGs 



Cluster Name 


Redshift 


X-ray luminosity 

(10" ergs ^} 


Cluster Mass 

(11)^ Hi .,1 


BCG (total) 


J-K 


Stellar Mass 

(10^-- M .,) 


XLSS J022303.0-043622 (J0223) 
XMMU J2235.3-2557 (J2235) 
XMMXCS J2215.9-1738 (J2215) 
RX J0848.9 + 4452 (J0849) 
RDCS J1252.9 -2927 (J1252) 


1.22 
1.39 
1.46 
1.26 
1.24 


1 l+°-l 

O Q+0.9 

•^••^-0.5 

6.6tl;l 


1.0 ±0.4 

3.1 ±0.7 
1.8 ±0.4 
1.8 ±0.4 
2.6 ±0.6 


17.72 ± 0.01 
17.34 ± 0.01 
18.72 ± 0.01 
17.00 ± 0.02 
17.36 ± 0.03 


1.82 ±0.01 
1.87 ±0.02 

1.83 ±0.02 
1.86 ± 0.03 
1.83 ±0.01 


0.61 ±0.08 
1.26 ±0.14 
0.39 ± 0.05 
1.30 ±0.15 
0.89 ±0.11 



The cluster X-ray luminosities are bolometric estimates taken from the literature and the cluster 
masses are M200 values (Supplementary Information). The errors on the cluster masses are based 
on the X-ray luminosity errors and the intrinsic uncertainty in the scahng relations. The J and Ks 
observations of J0223, J2235 and J2215 (Supplementary Fig. 1) were taken with the 8.2-m Subaru 
telescope and reach a 5a (s.d.) limiting magnitude J ~ 23.7 and Kg ~ 22.8 (23.1 in the case of 
J0223). The photometry for our data was calibrated using standard stars taken on the night in the 
Vega system. For comparison with previous observations we find that our J0223 BCG total Kg- 
band magnitude (Kg =17.72ib0.01) is in excellent agreement with the total magnitude from the 
literature {Kg =17.76±0.04, assuming a i^^ -band conversion from AB to Vega system of -1.86). 
The photometry for J1252 and J0849 were sourced from the literature ^"^'^^ and for these galaxies 
the total Kg magnitudes and J — Kg colours were measured in similar aperture sizes. All data 
have been analysed in an identical manner for direct comparison (see Supplementary Information). 
The errors on the stellar masses include all photometric errors and the uncertainty in the cahbration 
with the semi-analytic model ^. All errors are Icr (s.d.). For each cluster we identified the brightest 
galaxy from the Kg-hmd magnitudes of all galaxies within 500 kpc of the cluster X-ray centroid 
because for approximately 95% of clusters the BCG lies within this radius All identified BCGs 
have optical spectra confirming their cluster membership 
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Early assembly of the most massive galaxies 
Supplementary Information 



We provide further details of the observations and data reduction; carry out a com- 
parison of the cluster masses used in the Millennium Simulation with those of our 
high redshift sample; and we also discuss the BCG mass determinations in greater 
detail and how these are affected by uncertainties in the stellar formation history of 
the galaxies. 

Observations and Data Reduction 

The MOIRCS instrument on the Subaru Telescope provides imaging and low-resolution 
spectroscopy over a total field-of-view of 4' x 7' with a pixel scale of 0.117" per pixel. This 
is achieved by dividing the Cassegrain focal plane and then re-focussing the light through 
identical optics onto two HAWAn-2 2048 x 2048 CCDs, each covering 4' x 3.5'. Obser- 
vations were taken in photometric conditions in 0.5" seeing on the nights of August 8th 
and 9th 2007, with the clusters centred on Detector 2. A circular 1 1 -point dither pattern of 
radius 25" was used for both bands to ensure good sky subtraction. Integration times were 
25 mins at J and 21 mins at Kg (37 mins for J0223). These exposures reach a ba (s.d.) 
limiting magnitude of J ~ 23.7 and Kg ~ 22.8 (23.1 in the case of J0223). Supplementary 
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Fig. 1 shows combined J and Kg images of J0223, Jll^S, and J2215. 

The data were reduced using the external IRAF package MCSRED. The data were 
flat fielded, sky subtracted, corrected for distortion caused by the camera optical design, 
and registered to a common pixel coordinate system. The final reduced images on which 
we performed the photometry were made by taking the 3cr (s.d.) clipped mean of the 
dither frames. The BCG photometry was extracted in an identical manner to the compar- 
ison sample using SExtractor (version 2.5) 'Best' magnitude, which is found to be within 
0. 1 magnitude of the total. Measuring total magnitude is a significant improvement over 
the fixed aperture photometry of previous studies and constitutes the optimum compar- 
ison with semi-analytic models, which provide no information on the spatial distribution 
of light from merging haloes. Furthermore in the densely populated regions such as clus- 
ter cores the 'Best' magnitude initially reverts to the isophotal magnitude (MAG ISO), 
enabling light from close neighbours to be excluded, and then incorporates an empirical 
aperture correction. To calculate the colours of the BCGs we run SExtractor in dual mode 
so that the Kg band detections extract the J band catalogue with identical positions and 
apertures which ensures accurate colour determination. 
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Cluster Masses 

Various authors have identified a weak correlation between BCG mass and their host 
cluster mass 9.20,31,32 y^,^^^^ does not change significantly with redshift out to 2; ~ 0.8. If 
cluster masses have grown by a factor 2 or 3 since 2; ~ 1, as some authors suggest then 
their BCGs would be expected to grow by 20 — 30% over the same time interval. Therefore 
in order to compare BCG masses in a meaningful way it is necessary that our cluster 
sample is well matched to the masses of simulated clusters in the Millennium Simulation 
with which we are comparing. The clusters of interest from the simulation are the 125 
most massive systems at the two redshifts z — 1.08 and z — 1.5, selected for comparison 
with observations Halo masses M200 are measured at a radius (-R200) inside which the 
average mass density is 200 times the critical density of the Universe. 

In order to compare our data with simulations we use the bolometric X-ray lumi- 
nosity of each system (Lx) ^i i^- 34.35 jj^^g ^^j^ \x%q,6, to determine cluster mass from 
power-law scaling relations For J2235 we use the published X-ray luminosity in 
the energy range 0.5 — 2.0 keV and convert to bolometric luminosity assuming a tempera- 
ture of 6 keV. The mass-observable scaling relations are reasonably well calibrated at low 
redshift but are far less well measured beyond z ~ 0.5, where evolution becomes impor- 
tant. However, it has been shown 3^-39 jj^^j self-similar evolution provides a reasonable 
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description of the dynamical state of clusters and furt;hennore that simple luminosity is a 
reliable mass proxy, with an intrinsic scatter um — 21% (s.d.) out to 2; ~ 1. We adopt the 
self-similar scaling relation between Lx and M500 (where M500 is defined in an analogous 
way to M200) calculated for 115 clusters in the range 0.1 < 2; < 1.3 (ref. 38) and given by 



where C = 5.6 ± 0.3 x 10^^ erg s-\ a = 7/3, (3 = 1.96 ± 0.10 and E{z) describes 
the evolution of the Hubble parameter; 



Total X-ray luminosities for our clusters were measured using apertures of at least 50 
arcsec radius, corresponding to about 500 kpc atz — 1, which is close to i?5oo for massive 
clusters. From the Millennium Gas simulation the ratio of Lx at R200 and R^qq is only 
1.03, which adds weight to the assumption that Lx measured at i?5oo captures nearly all 
the X-ray emission. Finally, we calculate the M200 values using the conversion from M500 
(ref. 41). In the absence of cool core clusters at high redshifts there is no need to exclude 
the core emission from these estimates. The Lx values and final M200 mass estimates for 
our clusters are shown in Table 1 of the main text. 

To demonstrate stability we note that using cluster temperature, where available, as 




E{z) = [nmil + Zf + {l-nm- ^A)il + zf + J^a] 



1/2 
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a proxy for mass , gives virtually identical results. Independent mass measurements are 
available for some clusters which provide an important check on the self-similar scaling 
assumption. From a dynamical analysis of J2215 (ref. 18), the cluster velocity dispersion 

(Ty = 570 ± 190 km s~^ and i?2oo = 0.63 ± 0.15 Mpc. Adopting the relationship 

3 

gives M200 — l-4loi ^ 10^^ -^0' which compares very well with the value in Table 1 
of 1.8 (±0.4) X 10^^ Mq. From a weak-lensing analysis of the Lynx-East cluster J0849 
M200 = 2.0 (±0.6) X lO^^M©, compared to our value of 1.8 (±0.4) x lO^^M©. Using 
the X-ray surface brightness of J 1252 measured with Chandra and XMM -Newton gives 
a mass M500 — 1-9 (±0.3) x 10^'^ Mq, which compares favourably with our estimate 
measured at R500 of 1.9 (±0.6) x 10^^ Mq, corresponding to M200 = 2.6 (±0.6) x 10^^ Mq. 
This mass also agrees with the estimate within R500 from the weak-lensing analysis of the 
cluster Finally, an extensive optical spectroscopic survey of J 1252 (ref. 14) reveals a 
mass inside R500 of 1.6 — 2.3 x 10^^ Mq, again consistent with our estimated value. We 
therefore conclude that our mass estimates are reasonably reliable and consistent between 
independent methods despite the high redshift of the sample. 

Crucially these cluster masses are comparable to the massive haloes seen in the Mil- 
lennium Simulation. A histogram of halo and cluster masses is shown in Supplementary 
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Fig. 3. The simulated cluster samples at z — 1.08 and z — 1.5 have lower mass limits at 
these redshifts of 1.5 x 10^^ M© and 1.0 x 10^^ M© respectively. The average mass of the 
combined high redshift sample of 250 haloes is 2.3 (±1-1) x 10^^ M© (s.d.), compared to 
the average mass for our sample of 2.1 (±0.8) x 10^^ M© (s.d.). Although the predicted 
mass function is relatively steep and the halo numbers fall rapidly with increasing mass \ 
there are still 32 clusters at z — 1.08 and a further 6atz — 1.5 with a mass larger than our 
most massive cluster J2235. 

Stability of BCG Mass Estimates 

Most studies of the K band Hubble diagram for BCGs have assumed that the near- 
IR light is a direct proxy for stellar mass 9,20,45-47^ sometimes using colours to confirm 
the presence of an old stellar population. However, estimates of the stellar masses of 
galaxies are known to depend on the underlying stellar evolution model used, leading to 
degeneracy and systematic uncertainties. In particular some stellar evolution models '^^''^^ 
are known to produce younger ages and smaller stellar masses compared to others. In order 
to investigate these systematic effects in our BCG mass estimates, to provide plausible 
errors and to identify the major assumptions on which our results depend, we study the 
effect of using a representative set of different stellar models covering an appropriate range 
of values for the important parameters. 
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We concentrate on three simple-stellar population (SSP) synthesis models: Bruzual 
and Chariot (hereafter BC), which is based on the Padova stellar models and Geneva 
spectral libraries; Maraston"^^ (hereafter MAR) and BaSTI The last two stellar pop- 
ulation models are based on independent spectral libraries and stellar tracks and generally 
incorporate a larger range of parameters. However, the significant addition in the MAR 
and BaSTI models, compared to BC, is that they both implement improved modeling of the 
thermally pulsating phase of extremely bright asymptotic giant branch (AGB) stars. It is 
important for our purposes to include a model with a realistic prescription for AGB evolu- 
tion because their contribution to the near-infrared light can be significant, even though the 
stars are relatively few in number, giving rise to significantly smaller mass-to-light (M/ L) 
ratios and therefore smaller stellar masses for a given Kg luminosity. This AGB evolu- 
tionary phase is difficult to model in detail and can be a source of discrepancy between 
different population synthesis models, hence the inclusion of two independent predictions 
here. 

We consider the SSP models with a range of metallicities Z — 0.4 ^0-2. 5 Zq, where 
Z is the total mass fraction of elements heavier than helium measured in solar units, and 
four formation redshifts Zf — 2, 2.5, 3, 5. For BC we use a Chabrier stellar initial mass 
function (IMF), while for those of MAR and BaSTI we use a Kroupa IMF. These IMFs 
are similar and both account for the flatter slope of low mass (< 1 Mq) stars as observed in 
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star counts. The major difference between Chabrier and Kroupa is simply the parametric 
fit to the low mass slope of the IMF. 

We extract the appropriate M/ L ratios from the SSP codes for the above ages and 
metallicities, which gives us the variation in stellar mass for a given observed flux in 
the Kg band and therefore a realistic handle on the uncertainty of the BCG masses. We 
normalise the M/ L ratios to corresponding values ^ at z — and then compare the M/ L 
estimates using the appropriate age for the mean redshift of our sample (z — 1.3). In 
Supplementary Fig. 4 we plot the BCG stellar mass at z = 1.3 against metallicity and 
formation redshift by representing this quantity as a surface for each of the different stellar 
population codes and in Supplementary Table 1 we list each individual stellar mass. The 
range in stellar mass predicted by these models is about 0.3 dex. In general we find that 
the older populations have a higher mass-to-light ratio, hence higher stellar mass, than 
those with a more recent formation epoch due to the short-lived nature of massive bright 
stars. The masses derived from the EC and BaSTI models are found to be systematically 
higher at all metallicities and ages compared to those derived from MAR by ~ 0.1 dex. 
These results are reasonably consistent with those found by other authors ^9,53,54 
choice of stellar population model can therefore have a significant affect on the derived 
BCG masses, meaning the significance of our result can vary between: 0.32% — 1.11% 
for BC; 1.0% - 9.8% for MAR and 0.68% - 2.5% for BaSTI (Supplementary Table 1). 
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If we restrict the metallicity to Z > 1.0 Zq the significances remain greater than: 0.82%, 
4.42% and 2.53% respectively for Zf — 2.0. These significance Umits improve further if 
we let metallicity vary over its full range and instead restrict the formation epoch of these 
systems to Zf > 3, with values of 0.58%, 2.25% and 1.48% respectively. 

Turning to the observations, estimates of metallicities for BCGs in the literature 
are consistent with solar or super-solar abundancies Furthermore, in addition to the 
J — Kg colour in Fig. 2 of the main text, evidence that the stars in BCGs form early 
is supported by other investigations using the spectrophotometric properties of passively 
evolving BCGs and red galaxies at the centres of clusters. For example, a best fit model 
for the J0223 BCG from spectral template fitting using SSP models gives 2;/ ~ 3 (ref. 1 1); 
while a study of the CMR in the Lynx cluster (J0849) at 2; = 1.26 implies a mean 
stellar age of 3.2 Gyr, corresponding to Zf > 3.7, with the stellar content of the bright 
elliptical galaxies in place and formed hy z ^ 3. Similar conclusions are reached for 
the early-type galaxies in J1252 from two studies of the colour-magnitude relation, which 
both suggest Zf — 2.7 — 3.6 with subsequent passive evolution I'^ ss-sg conclude that 
since observations suggest Z > 1.0 Zq and Zf > 3.0, our central result of negligible 
stellar mass evolution is not seriously compromised by the age and metallicity caveats and 
remains significant for the most plausible stellar evolution histories of these galaxies. 
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In order to further investigate systematic effects, some authors ' have attempted 
to caUbrate stellar mass from the SSP models with independent estimates of gravitational 
mass from velocity dispersion measurements. This has uncovered potentially serious sys- 
tematic offsets (0.3 — 0.4 dex) between the two mass estimators for galaxies at 2; ~ 1. 
Such offsets are currently difficult to interpret as there is evidence that they may equally 
well be produced by biases in dynamical mass measurements or evolution in the dynam- 
ical structure of galaxies, as by biases in the photometric stellar mass determinations 
However, we note that the hybrid model shown in Fig. 2 of the main text (equivalent to 
an exponentially decaying model with r = 0.93 Gyr) is very close to the r = 0.97 Gyr 
model used in one of these comparisons which shows negligible offset (0.03 ± 0.06 dex) 
between the stellar and kinematic mass estimates. 
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Supplementary Figure 1: Infrared images of the three high redshift clusters observed 
with the Subaru telescope at the summit of Mauna Kea. Data were taken in the J and Kg 
wavebands using the MOIRCS (Multi-Object Infrared Camera and Spectrograph) instru- 
ment at the Cassegrain focus. The images show the 1.5' x 1.5' regions around the clusters, 
a, XLSS J022303-043622. b, XMMXCS J2235.3-2557. c, XMMXCS J2215.9-1738. 
The X-ray contours overlaid in green are taken from the XMM-Newton XCS pipeline, 
smoothed with a Gaussian kernel and normalised to the peak of the cluster emission. At 
these redshifts 1.5' corresponds approximately to 0.75 Mpc. 
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Supplementary Figure 2: The BCG Ks band Hubble diagram with the five high redshift 
BCGs shown in red. The photometric errors are ~0.01 mag as shown in Table 1 of the 
main text. The lines represent predictions of the observed Kg magnitude based purely on 
stellar population models with: no evolution (solid), passive evolution with 2;/ =2 (dotted), 
the hybrid passive evolution model with 50% of stars formed by 2; = 5 and 80% formed 
by z = 3 (dashed), equivalent to an exponential decay with r = 0.93. The black points 
indicate the comparison sample of 81 low redshift BCGs with masses in the same range as 
our high redshift sample. 
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Supplementary Figure 3: A histogram of the 250 halo M200 values from the semi- 
analytic model ^ ^/ = 1-08 and Zf = 1.5 (black) and the corresponding values for 
our observed sample (red). The average mass of the simulated and real cluster samples is 
2.3 (±1.1) X 10^4 Mq (s.d.) and 2.1 (±0.8) x IO^^Mq (s.d.) respectively. This demon- 
strates that our cluster masses are well matched to those in the models allowing for direct 
comparison of the BCGs. 
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Supplementary Figure 4: The three plots display the stellar mass surface obtained when 
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the metallicity and formation redshift are varied over the ranges Z — OAZq — 2.5Zq 
and Zf — 2, 2.5, 3, 5 respectively. The models considered are BC, MAR and BaSTI (up- 
per, middle, lower) as described in the Supplementary Information text. The grey shaded 
surface is the mean mass of our BCG sample while the black mesh surface represents a 
constant stellar mass of 5.55 x IO^^Mq corresponding to a significance for a one-tailed P 
of 0.05 (d.f . = 4) in the offset from the predicted BCG average mass. The masses derived 
from the BC and BaSTI models are found to be well in excess of this limit for all combi- 
nations of metallicity and age considered here. For the MAR model there is shown to be 
a steep drop in mass-to-light ratio for young stellar populations meaning this significance 
(0.05) is not achieved if both the metallicity Z < IZq and stellar population formation 
epoch Zf < 2.5 (see Supplementary Table 1). 
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Zf = 25 


Zf = 2 


BC 










I.OZ0 


0.97(0.54) 


0.96(0.58) 


0.93(0.63) 


0.91(0.71) 


O.4Z0 


1.10(0.32) 


0.96(0.58) 


0.93(0.63) 


0.82(1.11) 


2.5Z0 


0.97(0.54) 


0.98(0.52) 


0.94(0.62) 


0.88(0.82) 


MAR 










I.OZ0 


0.75(1.61) 


0.70(2.13) 


0.63(3.23) 


0.62(3.43) 


O.5Z0 


0.84(1.00) 


0.70(2.13) 


0.53(6.12) 


0.46(9.82) 


2.2Z0 


0.76(1.52) 


0.69(2.25) 


0.59(4.14) 


0.58(4.42) 


BaSTI 










I.OZ0 


0.92(0.68) 


0.80(1.48) 


0.77(1.44) 


0.70(2.13) 


O.5Z0 


0.91(0.71) 


0.84(1.00) 


0.82(1.11) 


0.76(1.52) 


2.OZ0 


0.86(0.91) 


0.80(1.23) 


0.76(1.52) 


0.67(2.53) 



Supplementary Table 1: The matrix of BCG steUar mass values (in units 10 M0) for our BCG sam- 
ple assuming z = 1.3 derived for a set of formation redshifts (zf) and metallicities for the three stel- 
lar population codes BC, MAR and BaSTI, as described in the Supplementary Information text. The 
numbers in brackets are the significance percentages (d.f . = 4) of the masses compared to the average 
(0.192 (±0.038) X 10^'^ Mq) of the 250 simulated BCGs. The stellar code used in the semi-analytic model^ 
is based on BC and almost identical to that shown in Fig. 2 of the main text; this gives a stellar mass of 
0.899 (±0.082) X 1O^^M0. Comparing the J - Kg colours of these models with the data we find that for 
Zf <2 with either solar or sub-solar metallicity the BC models are too blue at the 3cr level compared to the 
average colour of our sample (J — Ks = 1.84 ± 0.05, s.d.). Otherwise the colours are reasonably consistent 
providing no additional significant constraint. -3-3 
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